First-principles theory of magnetically induced ferroelectricity in TbMn0 3 
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We study the polarization induced via spin-orbit interaction by a magnetic cycloidal order in 
orthorhombic TbMnOs using first-principle methods. The case of magnetic spiral lying in the b- 
c plane is analyzed, in which the pure electronic contribution to the polarization is shown to be 
small. We focus our attention on the lattice-mediated contribution, and study it's dependence on 
the Coulomb interaction parameter U in the LDA+(7 method and on the wave-vector of the spin 
spiral. The role of the spin-orbit interaction on different sites is also analyzed. 

PACS numbers: 75.80.+q,77.80.-c 



INTRODUCTION 

Magnetically induced ferroelectricity provides a route 
to materials with a large magneto-electric (ME) cou- 
pling. The appearance of the ferroelectric order simul- 
taneously with the onset of a particular magnetic or- 
der suggests that the magnetic and electrical proper- 
ties should be strongly interconnected. Therefore, one 
might expect to observe a strong dependence of the elec- 
tric polarization on external magnetic fields in such ma- 
terials. Indeed, a pronounced interplay between fer- 
roelectricity and magnetism was observed experimen- 
tally [EH SI in (Gd,Tb,Dy)Mn0 3 , (Tb,Dy)Mn 2 5 , 
N13V2O8, etc. In these materials, ferroelectricity is in- 
duced by a specific (e.g., cycloidal) magnetic ordering. 
Since polarization is a polar vector quantity, it is neces- 
sary to break the inversion symmetry in order for the fer- 
roelectric phase to appear. In magnetically induced fer- 
roelectrics, the inversion symmetry can be broken by the 
magnetic order alone, or in some cases, by the magnetic 
order in combination with the crystal structure (e.g., in 
CaaCoMnOe @), where separately both the lattice and 
magnetic structures have inversion symmetry. The spin- 
orbit (SO) coupling is essential in order to allow the bro- 
ken symmetry to be communicated from the magnetic to 
the charge and lattice degrees of freedom. 

In our work, we focus our attention on TbMnOs, in 
which the appearance of polarization is attributed to the 
cycloidal order of magnetic moments on the Mn 3+ sites. 
TbMnC>3 is an orthorhombically distorted perovskite ma- 
terial whose lattice has Pbnm symmetry with 20 atoms 
per unit cell. The magnetic properties of TbMnOs are 
mainly due to the Mn 3+ ions. The competition between 
the nearest-neighbor and next-nearest-neighbor spin in- 
teractions leads to a rich phase diagram. Experimentally 
it is known that below ~ 41 K the magnetic moments 
form a sinusoidally modulated antifcrromagnctic order. 
Below ~ 27 K, a phase transition occurs in which the 
magnetic order develops a cycloidal character (with spins 
in the b-c plane) with modulation wave- vector k s ~ 0.28 
along b, and a polarization simultaneously appears along 
c. Phenomcnological and microscopic theoretical mod- 



els ((J 0, H| can explain how the cycloidal magnetic order 
can drive the system to become ferroelectric, and one 
can predict the direction of the polarization (but not its 
sign) on rather general symmetry grounds. However, it 
is unclear whether one should expect the ionic displace- 
ments (phonons) to play the most significant role in cre- 
ating electric polarization, or whether a purely electronic 
mechanism could explain the observations. 

Our previously-reported initial findings [l8| and a re- 
lated study [lj| suggested that, at least in the case of 
the spin spiral lying in the b-c plane, the lattice contri- 
bution to the polarization dominates over the electronic 
contribution in TbMnOs. However, these works left sev- 
eral questions unanswered. Here, we first briefly review 
our earlier work, emphasizing a careful analysis of the 
lattice-mediated contribution to the polarization using 
first-principles calculations. Wc then discuss several ex- 
tensions of the work, including a study of the dependence 
of the results on the choice of U parameter, additional 
details concerning the site-specific spin-orbit interaction 
and its effects on the dynamical effective charges, and a 
more thorough and revealing investigation of the depen- 
dence of the polarization on the wave-vector of the spin 
spiral. 



COMPUTATIONAL DETAILS 

The electronic-structure calculations are based on the 
density functional theory (DFT) and are performed u sing 
the projector-augmented wave (PAW) method 0, Elf 
in a plane- wave basis set. The plane- wave energy cut- 
off is 500 cV. We use the VASP code package [9( for 
our calculations. We do not include / electrons in the 
Tb PAW potential. The local-density approximation for 
the exchange-correlation functional is used (Ceperley- 
Alder [13] with Vosko-Wilk-Nusair correlation [l3|). For 
a proper treatment of Mn d electrons, we use on- 
site Coulomb corrections implemented in a rotationally- 
invariant LDA+t/ formulation We consider two val- 
ues of U, 1 and 4 eV. The results of a systematic s tudy 
with U = 1 eV were reported previously in Ref. [l8j ; 



2 




FIG. 1: Sketch of the a x 36 x c/2 orthorhombic cell of 
TbMnOa (Pbnm space group) showing MnOe octahedral tilts 
(top) and the cycloidal magnetic structure on the Mn 3+ sites 
(bottom). 

some of these will be reproduced here for comparison. 
We consider the cycloidal magnetic order having a spin 
spiral lying in the b-c plane and wave-vector k s propa- 
gating in the fr-direction (see Fig. [I). The experimen- 
tal value of the wave vector is incommensurate with the 
lattice, k s ~ 0.28 [H. Since we use periodic boundary 
conditions, we study commensurate cycloidal spin struc- 
tures with wave-vectors fc s = 0,1/4,1/3,1/2,2/3 and 1 
using appropriate supercells. The structural relaxation 
was performed with k s = 1/3 (60 atoms per cell), which 
is close to the experimental fc s (Fig.[T]). A 3x1x2 fc-point 
sampling is used. The Berry-phase approach is used 
for the calculation of the electric polarization. 



RESULTS 
Dependence on U parameter 

A crystal structure optimization was performed for 
the 60-atom supercell using U = 1 eV and {7 = 4 eV. 
These calculations were carried out without SO interac- 
tion to obtain centrosymmetric reference structures, with 
respect to which Berry phase polarization will be calcu- 
lated in the subsequent analysis. The lattice parameters 
and Wyckoff coordinates of the relaxed structures are 
given in Table HI For both values of {7, the calculated 
structural parameters arc very close to the experimen- 
tal ones, with the lattice vectors being in slightly better 
agreement with experiment for the case of U = 4 eV. 

The calculation of the polarization in the presence of 
the SO interaction for these relaxed structures yields 
an estimate for the pure electronic contribution to the 
polarization. We find P cloc = 32^C/m 2 and P clec = 
— 14^C/m 2 for U — 1 eV and U = 4eV respectively 



TABLE I: Experimental (Ref. [T3]) and theoretical (17=1 eV is 
taken from Ref. [3]) structural parameters for orthorhombic 
TbMnOs. 

Exp. U = leV U = 4eV 



Lattice vectors 


a (A) 5.293 


5.195 


5.228 




b (A) 5.838 


5.758 


5.775 




c (A) 7.403 


7.308 


7.343 


Tb 4c(x y 1/4) 


X 


0.983 


0.979 


0.980 




11 


0.082 


0.084 


0.084 


Mn 4b(l/2 0) 










Ol 4c(:r y 1/4) 


X 


0.104 


0.107 


0.111 




11 


0.467 


0.469 


0.465 


02 8d(x y z) 


X 


0.704 


0.699 


0.700 




y 


0.326 


0.320 


0.323 




z 


0.051 


0.052 


0.053 



(with the direction of polarization parallel to the c axis, 
as was expected from symmetry). In both cases, the 
magnitude is much smaller than the observed value of 
~ 600/iC/m 2 [l|, indicating that it can be neglected for 
any reasonable value of U. 

To analyze thoroughly the role of the ionic displace- 
ments, we computed the Hellmann-Feynman forces which 
appeared on the ions after the SO interaction was turned 
on. Keeping in mind that only zone-center infra-red ac- 
tive modes can be responsible for the appearance of po- 
larization, we filtered out the forces that did not act 
on these modes. We found that the infra-red active 
modes that were left had dynamical dipoles along c as 
expected. There are eight such modes: three associated 
with Mn atoms in Wyckoff position 4b, one each for Tb 
and Ol atoms in Wyckoff position 4c, and three for the 
02 atoms in Wyckoff position 8d [2(| ■ To find the ionic 
displacements induced by the SO interaction, we used 
the force-constant matrix [l8[ calculated in the absence 
of SO using the 60-atom supercell. (We note however 
that one could have used the 20-atom unit cell for this 
calculation, with fc s = 0, since the force-constant ma- 
trix depends only weakly on the wave- vector.) Updating 
the ionic positions according to the calculated displace- 
ments and carrying out Berry-phase polarization calcu- 
lations for the resulting structures, we find the total po- 
larization induced by the SO interaction. The results for 
U = leV and U = 4eV are P tot = -467/iC/m 2 and 
ptot _ _218/iC/m 2 respectively. Knowing the ionic dis- 
placements, one can also find the contributions to the 
polarization from each mode by calculating the effective 
charges of the modes. The results of such calculations 
for both values of U arc given together with the forces in 
Table HI 

Note that the values for the forces in the second column 
of the table represent the forces acting on the 'modes' 
rather than ions themselves. Each mode is character- 
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TABLE II: Forces F (meV/A), effective charges Z (e), and 
contributions AP to the polarization (fiC/m 2 ) from IR-active 
modes. See text for the description of the conventions used to 
describe the modes. The values for AP* are calculated with 
the SO coupling turned off everywhere except on Mn sites. 



Wyckoff 




U = 


leV 




U = 


ieV 


mode 


F 


Z AP 


AP* 


F 


AP 


Tb 4c, z 


0.43 


7.47 


-94 


-73 


0.47 


-42 


Ol 4c, z 


2.26 


-6.82 


81 


69 


1.46 


15 


Mn 4b, x 


-7.04 


0.57 


-13 


-14 


-2.00 


-3 


Mn 4b, z 


-8.93 


7.46 


-248 


-232 


-3.86 


-94 


Mn 4b, y 


-2.94 


0.55 





-1 


-1.06 


2 


02 8d, x 


5.06 


0.08 


3 


2 


2.00 


2 


02 8d, y 


3.57 


0.23 


16 


9 


1.59 


6 


02 8d, z 


4.41 


-5.74 


-234 


-208 


1.37 


-87 








-489 


-448 


-201 



Site-specific spin-orbit interactions 

To better understand the role of the spin-orbit interac- 
tion, we studied how the behavior of the system changes 
depending on the strength of the SO interaction and on 
the presence of SO on various atomic sites. We performed 
calculations of the forces with the SO interaction turned 
off on all sites other than Tb, then Mn, then O. Us- 
ing the force-constant matrix, we estimated the lattice 
contributions to the polarization and found them to be 
— 11/zC/m 2 , — 447/xC/m 2 and — 8^C/m 2 , respectively. 
This result shows that the SO interaction on the Mn 
sites is responsible for almost all of the lattice-mediated 
contribution to the polarization (see the column for AP* 
in the Table HI]) . We also calculated the purely electronic 
contributions to the polarizations for these three cases. 
Interestingly, the electronic contribution comes almost 
entirely from the spin-orbit effect on the Tb sites. Sev- 
eral calculations of the forces with modified SO strength 
confirmed that they depend linearly on the SO strength. 



ized by a unit vector in the (20 x 3)-dimcnsional space. 
Therefore, to get the actual number for the force acting 
on a particular ion, one must multiply the corresponding 
value from the table by an appropriate unit vector. This 
procedure will give absolute values for the ionic forces 
that are two times smaller than the numbers in the table 
for the first five modes, and 2\/2 times smaller for the 
rest of the modes. 

The decomposition of the lattice-mediated polariza- 
tion into mode contributions is discussed in detail in Ref. 
[lit ■ Here we compare the results obtained with different 
Coulomb interaction parameters. The results in Table HTl 
for U = leV and U = 4eV may seem different, but in 
fact they are qualitatively similar. The forces may be 
viewed as vectors in an 8-dimensional space, and one can 
calculate the angle 9 between them to find cos 9 ~ 0.98. 
Thus, the direction of the forces is almost the same, 
although the magnitude is different. This means that 
the underlying physical mechanism of the magnetically 
induced polarization does not depend strongly on the 
choice of U, while the magnitude of the effect does change 
with U . 

In the subsequent sections we will consider only U = 
1 eV, which gives better agreement with the experimental 
polarization. Moreover, in this case the calculated band 
gap matches the experimental value of ~ 0.5 eV jl5||. In 
general, a better strategy for choosing the U value would 
be to calculate the exchange parameters and fit those, 
rather than the band gap, to the experimental data. In 
the present case, however, such an approach leads to a 
similar choice [19( of parameters. 



Dependence on wave-vector 

Microscopic theoretical models involving the displace- 
ments of ions 0, [H, [23[ show that the Dzyaloshinskii- 
Moriya (DM) interaction can induce ferroelectric dis- 
placements of ions. Usually, only the interaction be- 
tween the nearest-neighbor transition metals is consid- 
ered in such models. As a consequence, the polar- 
ization is expected to depend sinusoidally on the an- 
gle between the spins of the nearest-neighbor Mn sites, 
P oc e n>n i x (S„ x S„') [iHj and thus sinusoidally on 
k s . To study the dependence of the lattice contribution 
to the electric polarization on k s , we have carried calcu- 
lations of the SO-induced forces for a 40-atom supercell 
(fc s =l/2) and an 80-atom supercell (fc s =l/4). Note also 
that the same 60-atom supercell can be used to set up 
a spiral with the wave-vector fc s =2/3. In general, if the 
supercell consists of n primitive cells, one can construct 
spirals with wave-vectors m/n, where < m < n. We 
also used the primitive 20-atom cell for the calculations 
with k s = and k s = 1. In all these calculations we kept 
the same structural coordinates as for the 60-atom struc- 
ture to make sure that we only change one parameter (k s ) 
in this study. We calculated the forces, and again filtered 
those that were IR-active. We find that the pattern of the 
IR forces matches almost exactly that of the 60-atom cell, 
i.e., the directions of the eight-dimensional vectors almost 
coincide in all cases. Using the force-constant matrix cal- 
culated on the 60-atom supercell and the effective charges 
from Table HH we can estimate the polarization for the 
new wave- vectors. The result is shown in Fig. [2l Actual 
calculations were performed only for non-negative values 
of k s , as symmetry arguments show the that polarization 
must be an odd function of spin-spiral wave-vector. The 
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density functional theory within the LDA+U framework. 
We compare the results for two values of U (1 and 4 eV), 
and find that the mechanism of magnctoelectric coupling 
is quite independent of U. The dependence of the polar- 
ization on the spin-spiral wave period is studied in detail. 
We find that it deviates significantly from the sinusoidal 
dependence expected from simple models. The polariza- 
tion is almost linear in wave- vector for absolute values of 
k s up to 0.5. 

This work was supported by NSF Grant No. DMR- 
0549198. 



FIG. 2: Dependence of polarization (circles, scale at left) and 
total energy (squares, scale at right) on the spiral wave-vector 
k s . The triangle indicates the extrapolation of the polariza- 
tion to the experimental wave- vector of k B — 0.28. 



values fc s = and fc B = 1 correspond to a collincar spin 
arrangement, in which cases the polarization was zero 
as expected. One can see that the dependence of the 
polarization on fc s deviates significantly from a simple si- 
nusoidal form. Surprisingly, the polarization is almost 
linear in k s up to k s = 1/2; such a linear dependence is 
expected in the long- wavelength limit, but that is a rather 
stretched assumption for fc s up to 1/2. This result indi- 
cates that nearest-neighbor DM models oversimplify the 
mechanism of the polarization induction, and that taking 
the next-nearest-neighbor interactions into account may 
be important. 

The experimental fc s lies in the range where we can as- 
sume a linear dependence. The extrapolation to fc s =0.28 
yields a value for the lattice contribution of the polariza- 
tion of about — 410/LtC/m 2 . Fig. Q] also shows the depen- 
dence of the total energy per formula unit on the wave- 
vector. We remind the reader that in these calculations 
the structural parameters were kept fixed, with only the 
directions of the magnetic moments on Mn sites being 
changed. However, one can still see that the wave- vector 
at which the minimum occurs is close to the experimental 
value of fc a . 



SUMMARY 

We have analyzed the lattice contribution to the 
electric polarization in the cycloidal-spin compound 
TbMnOs, with the spin spiral lying in the b-c plane, using 
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